Reporte de avance

Preparación

Librerías

library(ncdf4)

Cargar datos

Predictores

## Flujo de radiación, onda larga
dlwrf <- nc_open("../../data/train/dlwrf_sfc_latlon_subset_19940101_20071231.nc")
## Flujo de radiación, onda corta
dswrf <- nc_open("../../data/train/dswrf_sfc_latlon_subset_19940101_20071231.nc") 

Datos de las estaciones

stations <- read.csv("../../data/station_info.csv")

Resultados de entrenamiento

energy_train <- read.csv("../../data/train.csv")

Exploración de los datos

Todos los predictores tienen el mismo formato:

3 variables:

• intTime
• intValidTime
• <Variable de interés>

Cada variable tiene 5 dimensiones:

• time
• lat
• lon
• ens
• fhour

print(attributes(dlwrf$var))

$names
[1] "intTime"                     "intValidTime"                "Downward_Long-Wave_Rad_Flux"

print(attributes(dswrf$var))

$names
[1] "intTime"                      "intValidTime"                
[3] "Downward_Short-Wave_Rad_Flux"

# Ver el formato
print(dlwrf)

File ../../data/train/dlwrf_sfc_latlon_subset_19940101_20071231.nc (NC_FORMAT_NETCDF4):

     3 variables (excluding dimension variables):
        int intTime[time]   (Contiguous storage)  
            long_name: time as an integer (YYYYMMDDHH)
        int intValidTime[fhour,time]   (Contiguous storage)  
            long_name: valid time as an integer (YYYYDDMMHH)
        float Downward_Long-Wave_Rad_Flux[lon,lat,fhour,ens,time]   (Chunking: [6,3,1,4,1705])  (Compression: shuffle,level 4)
            _FillValue: 9999
            units: W m-2
            long_name: Downward_Long-Wave_Rad_Flux_Average (Average for  Mixed Intervals) @ surface
            cell_methods: time: mean
            GRIB_param_discipline: Meteorological_products
            GRIB_param_category: Long-wave_Radiation
            GRIB_param_name: Downward_long_wave_rad_flux
            GRIB_generating_process_type: Forecast
            GRIB_param_id: 2
             GRIB_param_id: 0
             GRIB_param_id: 5
             GRIB_param_id: 192
            GRIB_product_definition_template: 8
            GRIB_product_definition_template_desc: Average, accumulation, extreme values or other statistically processed value at a horizontal level in a time interval
            GRIB_level_type: 1
            GRIB_level_type_name: surface
            GRIB_interval_stat_type: Average
            GRIB_VectorComponentFlag: easterlyNortherlyRelative

     5 dimensions:
        time  Size:5113 
            long_name: Time
            units: hours since 1800-01-01 00:00:00
            axis: T
        lat  Size:9 
            long_name: Latitude
            standard_name: latitude
            units: degrees_north
            actual_range: 31
             actual_range: 39
            axis: Y
        lon  Size:16 
            long_name: Longitude
            standard_name: longitude
            units: degrees_east
            actual_range: 254
             actual_range: 269
            axis: X
        ens  Size:11 
            long_name: ensemble
            standard_name: ensemble
            axis: E
            description: 0 is control, other values are perturbation numbers
        fhour  Size:5 
            long_name: forecast hour
            units: hours

    7 global attributes:
        Conventions: CF-1.0
        title: Subset of data from 2nd-generation multi-decadal ensemble reforecast generated from the NCEP Global Ensemble Forecast System, mimicking version operational at NCEP/EMC circa mid-2012.
        institution: NOAA Earth System Research Laboratory (ESRL)
        source: NCEP GFS v 9.01, T254L42.  Control initial conditions from CFSRR.  Perturbed initial conditions from ETR.  Model error simulated with STTP.
        references: http://www.esrl.noaa.gov/psd/forecasts/reforecast2/index.html
        history: Subset created 2013-01-15 19:20:02 UTC
        comment: Original dataset generated on DOE's supercomputers at Lawrence Berkeley Laboratory through ALCC/ASCR grant.

Datos de interés

Cada variable tiene exactamente las mismas cinco dimensiones. Del resumen impreso arriba podemos rescatar los siguientes valores importantes:

Dimensión	Tamaño	Valor Mínimo	Valor Máximo
Latitud	9	31	39
Longitud	16	254	269
Ensemble	11	1	11
fHour	5	12	24
Tiempo	5113	1700586	1823256

Latitud y Longitud

Representan valores en la zona geográfica de interés.

En esta área, solo toman los siguientes valores:

• Latitud: [31° - 39°]
• Longitud: [254° - 269°]

Mapa de las estaciones meteorológicas y los sitios de medición

En estos datos, la longitud está en grados positivos desde el meridiano cero, llegando hasta 360°. Por otro lado, las coordenadas de las estaciones meteorológicas tienen este dato en valores que van de 0° a 180°, ya sea al Este u Oeste.

head(stations)

Esto será tomado en cuenta al etiquetar la longitud en los datos de entrenamiento.

# Primero obtener las variables de interés
dlwrf.data <- ncvar_get(dlwrf, "Downward_Long-Wave_Rad_Flux")
dswrf.data <- ncvar_get(dswrf, "Downward_Short-Wave_Rad_Flux")

# Etiquetar las dimensiones Latitud y Longitud
# Latitud: 9 [Índice 2]
# Longitud: 16 [Índice 1]
dim(dlwrf.data)

[1]   16    9    5   11 5113

# Valores obtenidos del resumen
lat.inicial = 31
lat.final = 39
# Se resta 360 para convertir de grados positivos a la representación que tienen
# las coordenadas de las estaciones.
lon.inicial = 254 - 360
lon.final = 269 - 360

dimnames(dlwrf.data)[[2]] <- seq(from=lat.inicial, to=lat.final)
dimnames(dlwrf.data)[[1]] <- seq(from=lon.inicial, to=lon.final)

dimnames(dswrf.data)[[2]] <- seq(from=lat.inicial, to=lat.final)
dimnames(dswrf.data)[[1]] <- seq(from=lon.inicial, to=lon.final)

# Muestra de como queda etiquetado
print(dlwrf.data[,,1,1,1])

          31      32      33      34      35      36      37      38      39
-106 247.019 234.019 211.019 196.019 192.019 186.019 179.019 177.019 182.019
-105 241.019 233.019 217.019 203.019 197.019 187.019 181.019 190.019 190.019
-104 245.019 245.019 231.019 220.019 216.019 207.019 199.019 210.019 204.019
-103 256.019 247.019 232.019 229.019 227.019 220.019 214.019 218.019 213.019
-102 258.019 249.019 242.019 235.019 231.019 228.019 225.019 224.019 224.019
-101 269.019 255.019 251.019 243.019 238.019 234.019 234.019 232.019 232.019
-100 291.019 262.019 260.019 255.019 247.019 240.019 240.019 237.019 237.019
-99  327.019 280.019 266.019 261.019 251.019 245.019 244.019 240.019 240.019
-98  340.019 300.019 275.019 269.019 256.019 249.019 247.019 242.019 243.019
-97  349.019 330.019 290.019 277.019 268.019 259.019 249.019 243.019 246.019
-96  351.019 343.019 310.019 296.019 278.019 277.019 262.019 250.019 246.019
-95  337.019 337.019 321.019 314.019 292.019 279.019 270.019 261.019 255.019
-94  324.019 335.019 324.019 329.019 305.019 293.019 279.019 267.019 262.019
-93  336.019 336.019 331.019 330.019 314.019 310.019 302.019 282.019 266.019
-92  338.019 334.019 332.019 329.019 319.019 314.019 312.019 301.019 280.019
-91  323.019 328.019 326.019 322.019 319.019 316.019 311.019 310.019 297.019

Ensemble

De acuerdo con el glosario de la agencia que provee el conjunto de datos, un ensemble se define como:

Colección de modelos numéricos que muestran resultados ligeramente diferentes.

e1 <- dlwrf.data[,,,1,]
e2 <- dlwrf.data[,,,2,]
e3 <- dlwrf.data[,,,3,]
e4 <- dlwrf.data[,,,4,]
e5 <- dlwrf.data[,,,5,]

comparativo <- c(e1[1,1,1,1], e2[1,1,1,1], e3[1,1,1,1], e4[1,1,1,1], e5[1,1,1,1])
print(comparativo)

[1] 247.0190 246.0190 247.5191 248.0089 247.2623

El primer ensemble es el de control, y es el que será utilizado para este modelo.

dlwrf.data.e1 <- dlwrf.data[,,,1,]
dswrf.data.e1 <- dswrf.data[,,,1,]

dim(dlwrf.data.e1)

[1]   16    9    5 5113

Forecast Hour

Es la hora a la cual la estación meteorológica hace su pronóstico sobre la variable de interés. Esto sucede cinco veces en el día: Empieza a las 12, y se repite en intervalos de 3 horas, terminando a las 24.

dimnames(dlwrf.data.e1)[[3]] <- namedSeq("fhour_", 24, steps=3, start=12)
dimnames(dswrf.data.e1)[[3]] <- namedSeq("fhour_", 24, steps=3, start=12)

print(dlwrf.data.e1[1,1,,1])

fhour_12 fhour_15 fhour_18 fhour_21 fhour_24 
 247.019  239.000  242.000  262.019  263.019 

Tiempo

Es la fecha en la que se tomó la medición de la variable. Esta columna corresponde directamente con una de las respuestas en los datos de entrenamiento.

Sus valores son obtenidos de la variable time contenida en el conjunto de datos.

# Todas las variables manejan el mismo tiempo
time <- ncvar_get(dlwrf, "time")

# Diferente formato
head(time)

[1] 1700568 1700592 1700616 1700640 1700664 1700688

head(energy_train$Date)

[1] 19940101 19940102 19940103 19940104 19940105 19940106

Comparando la variable time con las fechas de los datos de entrenamiento, se observa que su formato no coincide. Cada una está representando su tiempo de forma distinta:

Variable	Formato
time	Horas desde 1800-01-01 00:00:00
energy_train$Date	Fecha formato YYYYMMDD

Por lo tanto es necesario convertir nuestro tiempo en horas, a fecha.

dates <- as.Date(as.POSIXct(time*3600, origin='1800-01-01 00:00:00', 'UTC'))
dates <- gsub('-', '', dates)
# Además, energy_train$Date tiene sus fechas como ints
dates <- strtoi(dates)

head(dates)

[1] 19940101 19940102 19940103 19940104 19940105 19940106

head(energy_train$Date)

[1] 19940101 19940102 19940103 19940104 19940105 19940106

Con estos valores, se pueden etiquetar nuestros predictores:

dimnames(dlwrf.data.e1)[[4]] <- dates
dimnames(dswrf.data.e1)[[4]] <- dates

head(dlwrf.data.e1[1,1,1,])

19940101 19940102 19940103 19940104 19940105 19940106 
247.0190 257.0190 231.8486 241.1632 242.1974 247.8309 

Preparar los subconjuntos de entrenamiento

Predictores

En este primer experimento se hará un modelo solo para el sitio de generación ACME.

Para identificar los datos que son útiles para este sitio particular se utilizan las coordenadas:

print(names(stations))

[1] "stid" "nlat" "elon" "elev"

acme <- stations[stations["stid"] == "ACME"]
print(acme)

[1] "ACME"       "34.80833"   " -98.02325" " 397"      

# Se redondea porque los archivos de los predictores tienen coordenadas enteras.
# Se tomará la medición más cercana.
acme.lat <- as.character(round(as.numeric(acme[2])))
acme.lon <- as.character(round(as.numeric(acme[3])))
print(c(acme.lat, acme.lon))

[1] "35"  "-98"

Se determinó que las coordenadas más cercanas a ACME son: (35, -98).

Con estos datos extraemos las mediciones correspondientes.

acme.dlwrf <- dlwrf.data.e1[acme.lon, acme.lat,,]
acme.dswrf <- dswrf.data.e1[acme.lon, acme.lat,,]

# Para simplificar este primer experimento, se decidió promediar las predicciones
# a lo largo del día
acme.dlwrf.means <- colMeans(acme.dlwrf)
acme.dswrf.means <- colMeans(acme.dswrf)
acme.final <- data.frame(
  Date = names(acme.dlwrf.means),
  dlwrf = acme.dlwrf.means,
  dswrf = acme.dswrf.means)
print(acme.final)

Respuesta esperada

Ya que se tienen los predictores, solo hace falta conocer cuánta energía se generó en ACME con esos valores de dlwrf y dswrf.

energy.acme <- energy[c("Date", "ACME")]
colnames(energy.acme) <- c("Date", "ACME")

# Hay que tener el mismo formato para combinarlo con nuestros predictores
energy.acme <- transform(energy.acme, Date=as.character(Date))
print(energy.acme)

Y se combina con los predictores en un solo Dataframe para poder crear modelos.

energy.acme <- merge(energy.acme, acme.final, by="Date")
print(energy.acme)

Creación de modelos

Regresión lineal simple

Usando un solo predictor: dlwrf

Observaciones:

• Se detecta un patrón no-lineal en la gráfica de Residuals vs Fitted
• La gráfica Q-Q indica que los residuales no están distribuidos de forma normal
• No se identifican puntos de apalancamiento
• La efectividad del modelo es pésima, lo que era de esperarse
• El predictor dlwrf parece ser significativo

lm.fit <- lm(ACME ~ dlwrf, data=energy.acme)
plot(lm.fit)

summary(lm.fit)

Call:
lm(formula = ACME ~ dlwrf, data = energy.acme)

Residuals:
      Min        1Q    Median        3Q       Max 
-18039113  -4582784    656225   5697096  16272029 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1371158     591168   2.319   0.0204 *  
dlwrf          47330       1777  26.639   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7375000 on 5111 degrees of freedom
Multiple R-squared:  0.1219,    Adjusted R-squared:  0.1217 
F-statistic: 709.6 on 1 and 5111 DF,  p-value: < 2.2e-16

Observaciones

• Casi se elimina el patrón no-lineal en la gráfica Residuals vs Fitted
• La efectividad del modelo mejora marginalmente

lm.fit2 <- lm(ACME ~ dlwrf + I(dlwrf^2), data=energy.acme)
plot(lm.fit2)

summary(lm.fit2)

Call:
lm(formula = ACME ~ dlwrf + I(dlwrf^2), data = energy.acme)

Residuals:
      Min        1Q    Median        3Q       Max 
-18877647  -4515620    436038   5530706  16824318 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.771e+07  3.452e+06   8.027 1.23e-15 ***
dlwrf       -1.215e+05  2.188e+04  -5.554 2.93e-08 ***
I(dlwrf^2)   2.618e+02  3.381e+01   7.743 1.16e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7333000 on 5110 degrees of freedom
Multiple R-squared:  0.1321,    Adjusted R-squared:  0.1318 
F-statistic: 388.9 on 2 and 5110 DF,  p-value: < 2.2e-16

Regresión lineal múltiple

Verificamos que no haya correlación entre los datos

plot(energy.acme$dlwrf, energy.acme$dswrf)

lmm.fit <- lm(ACME ~ dlwrf + dswrf, data=energy.acme)

plot(lmm.fit)

summary(lmm.fit)

Call:
lm(formula = ACME ~ dlwrf + dswrf, data = energy.acme)

Residuals:
      Min        1Q    Median        3Q       Max 
-18874982  -1458088   1125228   2163414  21383542 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  82952.2   301536.3   0.275    0.783    
dlwrf        -5292.1     1005.2  -5.265 1.46e-07 ***
dswrf        52945.4      438.8 120.663  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3759000 on 5110 degrees of freedom
Multiple R-squared:  0.7719,    Adjusted R-squared:  0.7718 
F-statistic:  8645 on 2 and 5110 DF,  p-value: < 2.2e-16

lmm.fit2 <- lm(ACME ~ poly(dlwrf, 2) + poly(dswrf, 2), data=energy.acme)

plot(lmm.fit2)

summary(lmm.fit2)

Call:
lm(formula = ACME ~ poly(dlwrf, 2) + poly(dswrf, 2), data = energy.acme)

Residuals:
      Min        1Q    Median        3Q       Max 
-18628158  -1475913    955625   2254224  20281315 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      16877462      51776 325.969   <2e-16 ***
poly(dlwrf, 2)1 -36602605    4386138  -8.345   <2e-16 ***
poly(dlwrf, 2)2  33081936    3738991   8.848   <2e-16 ***
poly(dswrf, 2)1 507944582    4166306 121.917   <2e-16 ***
poly(dswrf, 2)2  40333521    3982481  10.128   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3702000 on 5108 degrees of freedom
Multiple R-squared:  0.7788,    Adjusted R-squared:  0.7787 
F-statistic:  4497 on 4 and 5108 DF,  p-value: < 2.2e-16